Unveiling the lamellar structure of the human cornea over its full thickness using polarization-resolved SHG microscopy

A key property of the human cornea is to maintain its curvature and consequently its refraction capability despite daily changes in intraocular pressure. This is closely related to the multiscale structure of the corneal stroma, which consists of 1–3 µm-thick stacked lamellae made of thin collagen fibrils. Nevertheless, the distribution, size, and orientation of these lamellae along the depth of the cornea are poorly characterized up to now. In this study, we use second harmonic generation (SHG) microscopy to visualize the collagen distribution over the full depth of 10 intact and unstained human corneas (500–600 µm thick). We take advantage of the small coherence length in epi-detection to axially resolve the lamellae while maintaining the corneal physiological curvature. Moreover, as raw epi-detected SHG images are spatially homogenous because of the sub-wavelength size of stromal collagen fibrils, we use a polarimetric approach to measure the collagen orientation in every voxel. After a careful validation of this approach, we show that the collagen lamellae (i) are mostly oriented along the inferior–superior axis in the anterior stroma and along the nasal-temporal axis in the posterior stroma, with a gradual shift in between and (ii) exhibit more disorder in the anterior stroma. These results represent the first quantitative characterization of the lamellar structure of the human cornea continuously along its entire thickness with micrometric resolution. It also shows the unique potential of P-SHG microscopy for imaging of collagen distribution in thick dense tissues.


Theoretical background of P-SHG
We derive here the polarization-resolved SHG signal from fibrillar collagen in the theoretical framework of tensorial nonlinear optics 1 . SHG process is described by a second-order nonlinear susceptibility tensor that induces a nonlinear polarization ⃗ ⃗, at 2 : Here i,j,k refer to the Cartesian components of the fields and ⃗ ⃗, ⃗ ⃗ is the fundamental field. We assume as in previous reports that the susceptibility tensor exhibits a cylindrical symmetry in the fibril frame (xyz), where x the fibril direction, and that the Kleinman symmetry applies 2,3 . While more complex tensors may be used, these assumptions have been shown to well reproduce P-SHG experimental data, while they limit the number of independent tensor components to only 2: and . We set , called the fibril anisotropy parameter, which has been shown to provide the orientation  of the peptide bonds to the triple helical axis: Let's consider a linearly-polarized fundamental laser beam at angle θ to the axis X in the microscope frame (X,Y,Z): Eq. (2) is written in the plane wave approximation and neglects the axial electric field component along the microscope axis Z due to strong focusing. We have indeed shown previously that strong focusing does not affect the determination of the in-plane orientation φ of collagen fibrils 9 . It only affects the contrast of the polarimetric diagram and impedes accurate measurements of the ratio of the two susceptibility components in the microscope frame ⁄ , which is not used in this study. The SHG signal intensity is then calculated as:

S3
Here are the components of the susceptibility tensor in the microscope frame (X,Y,Z). One gets 7,10-14 : where is the intensity of the fundamental laser beam, C is a parameter merging geometrical and other parameters, and: The SHG signal thus depends on the incident polarization angle and on the orientation of the collagen fibrils in the microscope frame, that is the Euler angles to change from the microscope frame (X,Y,Z) to the collagen frame (x,y,z): the collagen angle  in the imaging plane XY and the angle  out of the imaging plane. Since the collagen lamellae in the corneal stroma are composed of well-aligned fibrils, we consider here that all the collagen fibrils within the focal volume have the same orientation, in contrast to more complex approaches that introduce sub-micrometer orientation disorder 7,8 . Note that P-SHG is not sensitive to the polarity of collagen fibrils, like usual linear polarimetry, so that the  and  angles are defined in smaller ranges than Euler angles: , ∈ 90°, 90° .

Numerical processing of P-SHG
We use a FFT analysis to extract the collagen orientation φ from experimental P-SHG data in a fast and efficient way. Equation (4) is thus expressed as a function of its Fourier components: where: FFT processing of P-SHG data then provides the coefficients , and defined as: Our data-processing workflow thus provides 4 types of information in every pixel:  The SHG signal averaged over all linear polarizations: 〈 〉 . This <SHG> intensity is similar to the usual SHG image acquired with circularly polarized excitation and same total acquisition duration.
 The in-plane orientation  of the collagen fibrils, which corresponds to the first minimum of the polarimetric diagram ( Fig. 1.c). It is determined as the weighted average of the angles and 15 :  A coefficient of determination R 2 that compares the experimental data and the curve obtained from the FFT parameters , and in every voxel (Fig. S3). Indeed, the FFT analysis provide angles even if the experimental data do not match the theoretical equation (4) at all. Notably, the cylindrical symmetry does not apply in voxels encompassing crossing fibrils from sequential lamellae with different orientations. In this case, equation (4) is not valid anymore and higher order Fourier components are present in equation (6). We thus calculate: Angles determined with R 2 <0.7 are eliminated in the quantitative analyses and depicted as black pixels in the images. R 2 >0.7 are used to code the brightness in orientation maps.
 The anisotropy parameter in the microscope frame, which corresponds to the square root of the ratio of the 2 minima: / Fig. 1c and is calculated as: This parameter varies with and the out-of-plane angle  and can be used to estimate . The transverse reconstructions nevertheless provide  in a more precise and easy way in this study.
This P-SHG processing workflow is implemented using Matlab (MathWorks Inc.).     Figure S4: computation of the distance of 2 distributions A and B. The distance is calculated as the summation of the absolute difference between the distribution values at each angle: , ∑ |ℎ ℎ | . It is then normalized to the maximum distance between two normalized distributions, ie divided by 2.